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ABSTRACT 


A finite element formulation for solving axisymmetric 
transient heat conduction and thermal stress problems is 
developed in this thesis. The governing equations of 
mmecoupled, linear, isotropic thermoelasticity are discretized 
using quadratic isoparametric elements. A FORTRAN IV progran, 
using double precision arithmetic, is presented. Compact 
storage techniques for banded symmetric matrices are used. 

Comparisons between exact and computer solutions demon- 
Strate close agreement for a number of test problems. De- 


tailed instructions for using the program are included. 
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LT. INTRODUCTION 


Thermal stresses have become increasingly important in 
engineering practice during recent years. In power genera- 
tion higher cycle temperatures and use of nuclear fission 
wre largely responsible for this trend. This thesis describes 
secomputer program for finding temperatures and stresses in 
bodies having axisymmetric geometry and loading. The govern- 
migeequatvons are those of isotropic, uncoupled, quasi-static, 
iamear thermoclasticity. They are discretized by using the 
finite element method. A FORTRAN IV computer program using 
double-precision arithmetic has been written to solve 


problems of the following kinds. 


A. TEMPERATURE PROBLEMS 

The transient temperature vector, evaluated ae the nodal 
points, may be obtained for an axisymmetric body with the 
combination of insulated, convection, or constant temperature 
boundary conditions. For the convection thermal boundary 
condition, however, we may have fluid flowing with entry 
temperature prescribed as a linear function of time (RAMP). 
The program can handle up to 15 different ramps, each having 
@-ditferent flow velocity, and with discontinuities between 


successive ramps. 


B. STRESS PROBLEMS 
The program will generate load vectors for pressure load- 


PPM eCeNtrtupal fToadine, and axial force. Provision 1s 
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mage for direct input of one additional load vector. Stresses 


may be found for any combination of these loadings. 


C. THERMAL STRESS PROBLEMS 

Thermal stresses may be found for as many as 20 different 
temperature vectors which may be output of the temperature 
meanlem: or direct input.,.-In short, in this part. any combina- 


eaon,of the temperature and stress problems may be used. 
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TI. FINITE ELEMENT FORMULATION OF HEAT CONDUCTION 
IN AXISYMMETRIC BODIES 
A. METHOD OF FORMULATION 
For bodies of revolution under axisymmetric loading the 
mathematical problems presented are two-dimensional. The 


governing equation for non-steady heat conduction is 


kyl =. och, (1) 


| 


Where k is the thermal conductivity, p the density, c the 
specific heat, T the temperature and V the gradient operator. 
The superior dot denotes a time derivative. 


Applying Galerkin's principle [1] gives 


(aN VET! dvs fipe Nee Ted, (2) 
ees V S 


memes the antegral is over the volume V of the conducting 
body and N; isda shape Lunetvon used in representing the 
Eemnperature distribution. Jf S is the surface of the body 
and n the outward normal to surface, then using Gauss' 


theorem we can write 


=. = a oT 
: V (N. kVT) dV = N; lea dS. (3) 
Since 
foie kvl) dV = § (VN-)° (kV1) dv 
V V : 
+f N-V- (kVT) dV, (4) 
V 1 


weacan combine Eqs. 2, 3 and 4 to get 
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pcN, T dv + f (VN, )* (kVT) av = s N. kee ds. (5) 


Each node of the solid region has a separate discretized 
linear equation calculated from Eq. 5 using the appropriate 
shape function N.. Thus each of the volume integrals on the 
meot Nand side of Eq. 5 yields a square coefficient matrix 
gn the assembled set of equations. Calculation of these 
mmemrces 1S a Standard process, , Details are given by 
Zienkiewicz [l]. 


The discretized set of equations takes the form 


la 


is weibes wh; (6)” 


II-< 


where there is a term by term correspondence with Eq. 5. 
aieeteal Symmetric matrices C and Y represent, respectively, 
the thermal capacitance and thermal admittance. The elements 
wemine vector I are nodal temperatures. The vector v, dis- 
cussed in the following section, depends upon the thermal 
boundary conditions. 

In the present development piecewise constant material 
properties have been assumed, i.e., k, p and c are constant 
Within each element, but may vary from element to element. 
Also two-dimensional isoparametric elements are used. Appli- 


cable equations are summarized in Appendix A. 


in this text a double underline denotes a rectangular 
matrix and single underline denotes a column vector. 
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B. THERMAL BOUNDARY CONDITIONS 
Thermal boundary conditions affect only those scalar 
equations of Eq. 6 which correspond to boundary nodes. 
feeordinely>s: thesvector v 1s sparse. Also, in a single 
problem it is common to have different thermal boundary con- 
ditions on individual portions of the boundary. In what 
feltows the subvectors of v (distinguished by individual 
Superscripts) which correspond to separate boundary condi- 
erons are treated individually. 
im Insulated 
ties ucloare that for ansulated boundary conditions 
the subvector vy) of the right hand side of Eq. 6 corres- 
ponding to this boundary condition is zero, Since =— = 
Zo CONVECTION 
the heat transfer mechanism occurs in the interface 
Grechne Solid and®’fluid. If the fluid temperature is 6 and 
ene Neat transfer coefficient is h, then equating heat con- 
ducted away from the surface to the efflux from the solid 


[2] gives 
~k ta) = h(T-6), (7) 


where the subscript S means that the derivative is evaluated 
ac the surface. 


For constant h: 


ar : a 
pee rmpiace A Aya f Ped Se (8) 
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In what follows the fluid temperature 8 is taken to be a 
specified function of position and time. For purposes of 
discretization, the fluid temperature is specified at a 
discrete number of fluid "nodes." If 8 represents the 

fluid temperature at fluid node j, then the fluid temperature 


along the boundary may be represented by 


6°= 5 N:0. (9) 


where the N; are one-dimensional forms of the shape functions 


fceetor the solid. Substitution in Eq. 8 gives 


oT « 
Ne eG oe=) Bo Va- (02-1) -) 5 10 
ENG k B y;3 (0j-7)) (10) 


where the summation extends over the surface nodes and the 
* 
€Gerficients ye aresg@iven by 


2 
Ya; =h 2 NUN, dss. Eee 


Mesenblaing the contributions from Eq. 10, the subvector 


v (2) for the convection boundary condition may be written 


vf2) = y%9, (12) 


ghe contributions -f Yig3j from veq.. VO are anecluded by. 


supmentime -Ehe matrax Y (see Eq. 13 below). 
3. Constant Temperature 
If 6 represents the constant temperature desired at 
the wetted surface, then we can use the convection boundary 


20 


condition and replace h by a big number (say 10°.).. Since h 


is very large, then for thermal equilibrium the temperature 
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T at the surface will be forced to equal 6. So the subvector 
v 63) for the portion corresponding to the constant temperature 
Poundary condition can be obtained from Eq. 12. 

Upon the application of these boundary conditions in 
sesingle problem, the right-hand side vector v will be conm- 
bined from the corresponding subvectors and the finite ele- 


ment discretized equation becomes 


+ 


eee TiN, (13) 


+ * 
Mere Y = Y + YY". 


SeeeaxACT TIME SOLUTION WITH SPATIAL DISCRETIZATION 
Pemcous@dcr only the solution of Eq. 135 for v = constant 
aut; a at time t = 0. Let Ee be a particular solution 


(steady state) with ie = 0 so that 


1 


i: = cae v (14) 
For the homogeneous equation 
: + 
C letyt, T= 0 (15) 


ine assumption T = w exp(=r’t), where w is a vector and A is 


eescalar, vyaelds the form 
a= Ay Cow GED) 


It is apparent that Eq. 16 defines an eigenvalue problem. 
Let A be the spectral matrix and W the modal matrix with 


normalization according to 


Oe a (17) 
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where I is the identity matrix of the same order as Now 


la 


let 


T = W exp(-At) b (18) 


meee ) 1S a constant vector. Substituting this in Eq. 15 
gives 
ieeheexp (sit) be 3G W Avexp(-Ar) \b. G9) 


memeeedq. 19 is satisfied for all b if 


weytwea, (19") 


gma this iS puaranteed to be satisfied since A and W are 
=wcetral and modal matrices for the eigenvalue problem of 
Eq. 16 with W normalized according to Eq. 17. 

Returmame tothe original problem (Eq. 13), the com- 
prete solution may be written as 

eye ee exp (2h) b) (20) 
where 
eee | (21) 
and 8 is a constant vector. Now 8 may be found (using 
Eq. 14) to be 

ee Wk (22) 


PPsrtBeuerne tis result into Eq. 20 and using the initial 


condition gives the result 


BesiWa: avchad? Wo vis (23) 


The general solution of Eq. 13 may thus be written 


as 
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T = W(L- exp(-At)) Av’) wiv + W exp(-At)W ta. (24) 


For purposes of the present program non-zero compo- 
memes Of vector v.are to be specified,as piecewise linear 
Benections of time. During each segment of time history of v 
Beanalytical solution of Eq. 13 is possible in the form of 
a particular solution plus a complementary solution such as 
Fa. 20.. At each node the corresponding time variation of 
Pemperacure will consist of a linear part contributed by the 
particular ee and a sum of n terms representing the 
Gomplementary part. Each of these n terms decays exponen- 
pee wttliva separate time »constant. In principle it -is.a 
straightforward process to find each particular solution and 
accompanying complementary solution. 

Contemplated problems may typically have from 10 to 
40 segments required to represent the piecewise linear 
wamtacion Of v. The number of body nodes n will be of the 
eraer of 100 or more. In view of the number of particular 
solutions required, each accompanied by an individual comple- 
mentary solution of the form given by Eq. 20, it was concluded 
Eat a numerical solution of Eq. 13 would be considerably 
more economical than an analytic one such as that given by 


Eq. 24. 


D. TIME INTEGRATION 

iM thrs Section the relative merits of the Runge-Kutta 
and trapezoidal methods of time integration are discussed. 
Since either of these methods will give an exact result if 


Ene Ysolucion is\’a linear function of time, investigation 
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Gs confined to performance on a single scalar equation 


ata cy > 0 (25) 


whose solution y = Des exp(-At) is of the same form as the 
components of the complementary solution (Eq. 20 ). 
1. Runge-Kutta Method 
A method introduced by Runge and subsequently elabo- 
rated by Heun and Kutta [3] is widely used for the numerical 
Betucion of first order ordinary differential equations. 
airs algorithm prescribes a sequence of calculations for 


determining the ordinate Y; ae tame. t- 


= 7. +At in terms 
nil rise] Ta e 


of Ys and values of y at intermediate and end points of the 
maeerval At. The fourth-order form, which requires four 
euaiuations of y, gives for Eq. 25 the result 

it = 1 - At + ary? (ary? + ar)" : (26) 
The right-hand side of Eq. 26 represents the first five terms 
Sree Taylor expansion of the’exact solution 
egg /V5 = exp(-AAt)) so we may conclude that the relative 
tor in cach time step 1s less than modulus of the next 
term: Gute Raine 

In addition to providing the apparent prospect for 
high precision indicated by this error bound, the Runge-Kutta 
method also permits changes of time increment during the 
integration process without requiring additional computa- 
tionally expensive matrix decompositions. The attractiveness 
of these two features dieeaeed a thorough exploration of the 


potential of this method for the present application. The 
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is confined to performance on a single scalar equation 
y thy = 0 (25) 


whose solution y = A exp(-At) 1s of the same form as the 
components of the complementary solution (Eq. 20 ). 
1. Runge-Kutta Method 
A method introduced by Runge and subsequently elabo- 
macea by Heun and Kutta [3] is widely used for the numerical 
Solution of first order ordinary differential equations. 
amas algorithm prescribes a sequence of calculations for 


cetermining the ordinate Y; ae tamer © 


= —.+AtT in terms 
evil ceil) al 


of ye and values of y at intermediate and end points of the 
muverval At. The fourth-order form, which requires four 


Exatuations of y, gives for Eq. 25 the result 
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iicerieht-hand side of Eq. 26 represents the first five terms 
Of the Taylor expansion of the exact solution 

(ig / 75 = exp(-AAt)) so we may conclude that the relative 
Etaonr in each time Step is less than modulus of the next 
term: (uke) / 50. 

In addition to providing the apparent prospect for 
naph precision indicated by this error bound, the Runge-Kutta 
method also permits changes of time increment during the 
integration process without requiring additional computa- 
tionally expensive matrix decompositions. The attractiveness 
of these two features dictated dsthorough exploration of.the 


potential of this method for the present application. The 


(ag 


- | -s war we 
a | 
0 Gpenated -F 

4 
roitéeupe weieoe signuie's xo: 81 


~ ey =4 


7 . i 7 - : As ~ 
Rij7 -« ‘i of* HNhS vag i» er (rk-)qxe Py * ¥ LIM 


oT) celtaloe yeersasmefgnoa sit do leF 


: ee 
si -sgnum 


. wih S 
5 Nel 
rss ‘2 ug ij UIT apres i Ley rd Leaubor snk hoai om = 
“1c 7h) haan vVighiw €f £27 soe bee nerell 
nope (aisns et htb yrentinw reba Peek Se " 
“feo. V1 ‘supee © 2edivoende alt: 


ead + vr Sot? te 2. .v ieee ae 


3" teow Jerre) +abeo-adriet oat 


ripes if <4 J i zoveg , 4 Ge 


Faun iwvean o@¢fs 
iaiopss tomiv?w eessege acne 


ra itigegrins) sivaee etehegee a 


disqualifying defect which emerged after studying a number 
of examples is readily appreciated from examination of 
Mmaose 1. For values of AAt less than 0.5 it is apparent 
mhat Runge-Kutta scheme affords acceptable engineering 
accuracy. However, when the method is applied to solution 
eeeeq. 15 we must deal with a number of A's equal to the 
mumpber Of nodes (see Eq. 20 ). This number may be greater 
enan 100 and the ratio of the largest A to the smallest may 
easily exceed 1000. Although the solution is dominated by 
meme contributions of the eigenvectors corresponding to the 
erarler XS; it 15 Clear that the solution will be unstable 
meeenec largest AAT exceeds about 2:7. Because an unaccept- 
aolvesmall At is required in typical problems, the Runge- 
Kutta method was rejected. 

Z@_ trapezoidal Method 


The trapezoidal method estimates ere from the formula 


4, NGG . . 
LG] \ La o OP (Y; os Venae: (27) 


substituting for ie and ys fromerq. 25 and rearranging 


+1 


gives 


ele Dee AAT 


ya 2 AMT ee) 


BeELecuexpansiton Of the Tipht-hand Side provides an error 


bound (per step): 
5 
WAT) 712. 
From the point of view of the size of AAt this method has no 


Seapiiaty tamic, but has slow attenuation with alteration 


22 


7 
= 


4) gaiyhare ¢ag%e 8 


=| wera : 


® ‘ , 3 


ui 9) 0 ihom ody wake tevewok 


* 9 aa) geek 708 
i Jeastasuse iz0Mia mse i he 


_ 


tatins- 2 sotw\ seeh 1306 ay 
| re £09 . pi See}  aebow 
‘ecote) G09 Pe @heeq edad 
ibe Gite hed 4th ode 5 


i104 GH te Groitdae 


vltures a9, tel) Geeks ee oes A 
<a Jive n= ee rat fumes 
ty 4 out (VEGF 
‘ 4 
ee aT} be U Aas 
; Kind Conny ¢ 
Ly t 4 
. 
via? ion mn { 
I 
‘i ; _ 
earls ol 70 


TABLE I 
ATTENUATION FACTOR CCMPARISON 


Fourth Order Runge-Kutta Algorithm 


At Yael; eee eit) 


(Runge-Kutta) (Exact) 

.0001 .9999 .9999 
S001, * .9990 .9990 
701 .9901 .9901 
<i .9048 .9048 
2 .8187 ae 187 
.§ 6068 6065 
1.0 3750 3679 
20 FOURS .1353 
aes .6484 0821 
3.0 125750 .0498 
4.0 5.0000 * 20183 
8.0 11025535 0003 
10.0 291.0000 0000 
20.0 5514.3333 0000 
50.0 240784.3333 0000 
100.0 4004901.0000 .0000 


ZS 


ny 


1 timer = 7 
NORISATOD AOTSAR WOT sia 


aisiqeglA stiod-sgead teahao dt 


oo — 


4 maar eT a! 
(6.2200 -epoua] : 


> teem ramenner— lie canat 


“Oo On 


40 O08 


—_— 


) ot7 Tee 


zy 


, 


noe 


mesien for large AAt. Table II shows this behavior. 

Since a wide usable range of AAt is essential and 
the stability of trapezoidal integration is guaranteed, 
this method is heees for the present program. 


Applying the trapezoidal algorithm to Eq. 13 yields 


i Soe eae (29) 
where 

n= ee 

Be ee 

oo 


and the eiterscrivts denote evaluation, at discrete time 
macervals At. If mis the order of the capacitance and 
Zecance Matrices, C and Y, then once a certain step size 
Beets chosen, it requires m>/3 Operations £O perrorm the 
meeded triangular decomposition of A. Thus, fone lease Mi, «a 
change of step size At becomes costly from the point of view 
of computer time. Accordingly, in the present program only 
One time Step size is used throughout each problem. 

Miso, 820m assuring SUtticLlently rapid attenuation of the 
components corresponding to the large AAt, the following 
correction is utilized. 

a tons Correction 

Irons proposed a scheme [4] for augmenting the atten- 
uation of the contributions of those eigenvectors for which 


AAt is large. Define 
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TABLE I] 
ATTENUATION FACTOR COMPARISON 


Trapezoidal Integration 


a Yael 7i erty 
(Trapezoidal) (Exact) 

0001 9999 9999 
001 9990 9990 
=O .9901 9901 
i! .9048 -9048 
a2 Housi2 . 8187 
oS : ood .7408 
9 -6000 6065 
1.0 so55 3679 
2.0 0000 okSS5 
25.5 = eta: wOiSi2a 
$<.0 =~ 2 0010 0498 
4.0 =s 9505 OLS 
8.0 -.6000 0003 
10.0 =. 6667 0000 
20.0 Son st ey -0000 
50.0 eeu 0000 
100.0 -.9608 0000 
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where Yi and Yi4 are obtained from Yk by trapezoidal 
mmcegration. 

In the program presented in Appendix B, Eq. 30 is used 
aiter every 10 steps of time integration. Table III shows 


the resulting modifications. 


E. ESTIMATION OF EXTREME EIGENVALUES 
Analytic results for one-dimensional heat conduction 
mive, £0r an eigenvalue, 

,- T° | (31) 
where d is the distance between points of extreme tempera- 
ture and zero temperature. 

fiewe use this to estimate the smalllest X in cylindrical 
coordinates, two modifications are recommended. 
ionssume that che point of zero temperature is in 

the fluid at a distance from the wall equal to k/h, 

Where hows the surface heat transfer coefficient. 

Zz. If there are two approximately orthogonal paths for 
heat flow from the (single) maximum temperature point, 


then replace wie in the above formula by 


= == 5 des (32) 


2 2 
d ote eee 


For estimating the largest i, the surface heat transfer 


eeerticient Nas no significant effect. We may continue to 


26 


7 
“THA OF ietehelGq) magores 
enh) Dey anpete ohm : 


vote bom 9x 


o'! ote, 2 HO 
vol 49: 1@ee9 a4 


Bc BY re ‘> ag 4 


‘ 7 _ bn 


oe 


6 


PTSeML ho Tae 


an kil ¥ taist 5 


Prk, Lit 


EFFECTS OF USING IRONS' CORRECTION 
AFTER 10 STEPS. OF INTEGRATION 


AAT eS ae Y10/Y0 y"10/%o 
EXact Trapezoidal Corrected 

0.00010 0.99900 0.99900 0.99900 
0.00020 0.99800 0.99800 0.99800 
0.00100 0.99005 0.99005 0.99005 
0.01000 0.90484 0.90484 0.90486 
0.10000 0.36788 0. 56757 0.36849 
0.20000 02.15554 0.13443 0.13579 
0.30000 0.04979 0.04866 0.04978 
0.50000 0.00674 0.00605 0.00645 
1.00000 0.00005 0.00002 — 0.00002 
2.00000 0.00000 0.00000 ; 0.00000 
4.00000 0.00000 0.00002 -0.00001 
8.00000 0.00000 0.00605 -0.00040 
10.00000 0.00000 0.01734 -0.00072 
20.00000 0.00000 0.13443 -0.00136 
50.00000 0.00000 0.44914 -0.00072 
100.00000 0.00000 0.67028 -0.00027 
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use the same formula for + » but now consider only the 
d 
smallest element and take 


A Seven eh Oresmawwves tt. Side 


min 4 4 


(32") 
d S lengin or largest: side 


max 4 


Mcomparison Of .eStamates based on. Eq. 31, 32, 32' with 
the exact solution of the eigenvalue problem has been carried 
Out for several examples. Based on these comparisons, it is 
believed that these estimates are sufficiently accurate for 
enmoosing a time step and estimating the time of occurrence 


of the extreme stresses. 
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III. ONE-DIMENSIONAL HEAT CONDUCTION 


For comparison of numerical (time) integration methods, 
studies of one-dimensional heat conduction were made. In 
Bias S€CLION. numerical results for trapezoidal time integra- 
faon using lrons' correction are compared with the exact 
eeansient temperature solution. 

Gonsider a tlatysiab of thackness L. with. conductivity k, 
Mmeoteyep, Specii1cerheat, c,. zero initial. temperature, one 
mace insulated and the other in contact with fluid at tem- 
perature Tr Cisse l)} .  thessumtace heat transfer is h. The 


Bxaece transient heat conduction solution is available [5] as 


oo Sinn, L Cosu,x “us wt 
T=T 1-20 —-—.-;——-; 'e (35) 
rg Ae u,L+Sinn, Lb Cosu,L 
Y 
~ Fluid 


“Temperature Ty 


Figure 1. Slab with One Face Insulated and Another 
in, Contact wath Fluad.. 
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where k/pc is the diffusivity and the H; are the solutions 


of the transcendental equation 


= ol 
Tan uL : ta 


(34) 

For the finite element comparisons we subdivide the 
distance L into m one-dimensional 3-noded elements and use 
the corresponding isoparametric shape functions. (The work- 
imenequations, shape functions and the element capacitance 
and admittance matrices are given in Appendix A.) 

In the course of this investigation separate computer 
programs were written to evaluate nodal transient tempera- 
tures using the following methods: 


ta) exact) transient temperature solution, Eq. 33; 


(b) exact time solution with spatial discretization 
Gsectron’ Tiss Eq. 24) 3 


(c) Runge-Kutta time integration (section II.D, 
Pert 2); i 


tay) trapezoidal”~taime integration, (section II-D, 
part b); 


(ee) trapezoidal time integration with Irons' 
cornectron (section it); "part c): 


For an initial step change of fluid temperature from 
zero to 1, transient temperatures have been found. For these 
comparisons the parameters (in consistent units) were taken 


EO be: 
Peo — 25, kk = 18, ¢ = 5 and h = § 


The distance L was subdivided into m = 3 elements. For the 
Present Plurpose, Comparison is confined to the exact solu- 


Elone(1tem (a)), and the finally adopted system (item (e)). 


30 


oy 


rotisioe why sre Pha ond. 


; 
iivibivge sw 2noelteqmes atte wf anh a: t 
tl cans (ea Qebed-2 (naodetomhh-sae - a. 3 a 
foitoue: ayode oiviemervagos? getive 
ineqgéa Jnomete of floes eno ita? squeal: Sem 


ciLiieqqs of nevig wee eooitiem sone 


ial 


‘iv? tye raga. pv lingiseavaeh etd eo Seyveos J 
| _ 97nnave O2 osITiw 19 
S2hotream 74 ivatte? ant rs 


{. > oyitdeys:nes Inelenety Bee 

-_ 

te ) lottnga wttiw oolSnfes gale Gee a 
‘ od ,J.2t nolises 


1OLI0TGoTeT GREP OF 90i-a 
(a orm 


vulgar) Quit: Lebioseunes 
'' otras?) ees taateos 


aH ori 1© “greasy <322 inissat ae ‘ 
feta 2 ‘.19peey tasiengyh) | & 

2 cf) eeeteqeteg ‘std  eped] 

on 72 (8 eo@ 935)" -9,,8 2 


al rT bohividdua 2cyu.J ssanie? 
jdt NOD eg Retea TEGHOR , 920Tiog If 
© Yiimuns? ody See, ((s) werry 


A) 


’ 


In Table IV, temperatures at the two faces (x = 0, x = 8) 
are compared for various times. The trapezoidal integration 
has been performed using the constant time increment unity 
and the Irons" correction is applied after every 10 incre- 


ments. 


It is believed that Table IV demonstrates that the 
momerical integration method gives adequate accuracy for 


engineering applications. 


TABLE TV 
COMPARISON OF ONE-DIMENSIONAL TRANSIENT TEMPERATURES 


Time=10|Time=20]Time=30] Time=50!Time=70 


Exact 


Trapezoidal 
with Irons' 


Exact 


Trapezoidal 
with Irons'Y 


“Tages correction is used after every 10 steps of trapezoidal 
integration. 
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IV. STRESS PROBLEM 


For bodies of revolution deformed symmetrically under 
Ersoymmetric loading, the stress distribution is two- 
Bemensional. Since deformation is symmetric about the axis 
eeerevolution, cylindrical coordinates (R,Z,6) are used. 
Beetollows that the Stress components are independent of the 
emueve 6 and all derivatives’ with respect. to 6 are zero. 


miso;the components of shearing stress Tp and tT vanish 


8 Z6 
on account of the symmetry. But since any radial displace- 
ment induces a strain Eg in the circumferential direction, 
PEesenon-Zero Component Of Straim and the three in-plane 
components (E75 Ep» Ypz)> complete the state of strain at a 
Peimt in any axisymmetric Situation. Hence the state of 


stress for an axisymmetric body under axisymmetric loading 


mo piven by 


(ofa 6 ps a ca (35) 


AY RZ 


Db akia 0 
In this chapter .the stiffness, .matrix.of.an axisymmetric 
body and the thermal, pressure, and centrifugal load vectors 
gre formulated and, finally, evaluation of stresses at a 
point is discussed. The treatment closely follows that of 
Zienkiewicz [1] and this reference should be consulted for 


further mMetails. 
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fe OF IFFNESS MATRIX 

The elements used are bodies of revolution (about the 
Seaxts). For analysis,it is sufficient to describe the 
Seess- section in the R,Z plane. In Fig. 2. such an element 


Sade the local §,n coordinates are shown. 


Figure 2. Quadrilateral Element Representation. 


If u and w are the displacement components at a point in the 
armections of R and Z respectively, then these displacement 
components may be defined in terms of the nodal displacements 


by the appropriate isoparametric shape functions as 
(36) 


where N., AeUMetsT ONO ec. dnd.t) 4) a5 the shape.function for 


element node i and u;,W;, are the nodal displacement compo- 


jes 
nents. The strain-displacement relations [6] can now be used 


to obtain the components of the strain vector. Thus 
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where B is the standard rectangular strain-displacement 
matrix of any finite element formulation, a function of 
eee tocal coordinates € and n, and 6 is the vector of nodal 
G@isplacement. (See Appendix A, part 3, where the applicable 
formulas and useful equations are summarized). 

Peetne elasticity matrix for an isotropic material is 


memenen the stress vector o ati a‘point is given by 
o=-De. (38) 


How.) by evaluation of the total strain energy in the element, 


tne element stiffness matrix can readily be obtained as 


Ko = s BT DB av, (39) 


where the integration extends over the volume of the 
element . 

iDeche present program the upper triangle of each ele- 
ment stiffness matrix is evaluated by numerical integration 
aeampe tour Gauss points within the range of &— and of n [1]. 
The element contribution is immediately placed in the system 


Stat inessimatrix, whicheis’stored:in banded form. 


B. THERMAL LOAD VECTOR 
If we denote T as the difference between local tempera- 
ture and reference temperature, then the thermal strain Ey 


is given as 


Cts UaT, (40) 
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where 


Bbes< in In doo>? 


and a is the coefficient of thermal expansion. 


The thermal load vector Ee is given by Zienkiewicz [1] 
as 
F=f BY ne” dy (41) 
—T = EG) a 


From.the point of view of the numerical evaluation it is 


e€ 


interesting to note, however, that the product D «> in 
Eq. 41 will reduce to 
pe enchant 
wu 6 P=2y, Dies ta) 


mere E is the modulus of elasticity and v is Poisson's 


meero. Thus 


eS  38@! L 
Bt ery J Eg! us (43) 
n 
where T= 2 NT; and m1 Isr the number of nodes in each 
i=l 


eLement. 

In the attached computer program in Appendix B the advan- 
Pape Or the simplicity of Eq. 42 has been utilized. Also, 
Sance B has Some zero components, in the process of multi- 
plication of a U, simply the addition of the appropriate 
non-zero components of each column of the Be has been per- 


formed. Finally, Eq. 43 has been integrated with four Gauss 


points. 
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Pee PRESSURE: LOAD: VECTOR 

Considerfa quadrilateral’ element as in’ Fig:'3.on the 
boundary of the axisymmetric cross-section where constant 
Meessure P is applied. The infinitesimal force dF due to 
the normal pressure acting on the inner infinitesimal cir- 


eumrerential surface dS is 
dpe="PdS-= 27 RP dk, (44) 


merc at is the infinitesimal length along the side of 


quadrilateral. 


Z,w 


Figure 3. Boundary Element Under Pressure. 


Let FR and Fo be the components: of the pressure force in the 


R and Z directions respectively, then 


dF, ="ak iGos “6, dF, ='-dE Sin @ ; (45) 
Where Sin 6 = ea and Cos 6 = ce . Therefore 
dF, = (27 RP)dZ , dF, = -(27 RP)dR. (46) 
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Since the total work W done by the normal force F is equal 


to the sum of the work done in R and Z directions, 


W 


Su dF, + Sw dF. 


ZTE i heueadze SR iw dR). | (47) 


New, by using the appropriate shape functions, each compo- 
nent of Eq. 47 may be defined in terms of the nodal values. 


Since 


i 
w= (6°) FE, (48) 


where ssh rs the *element pressure load «vector 


meetor CONntributed by node j 1S explicitly given as 


e al 
: FR 1 Lae 43 
ine = 2nP f (ENR,) No des. (49) 
ePane | Bed . -1 aN, q 

Jp Sy a Re 


where N. Hime 49 15 the appropriate shape function for 


node j. 


D. CENTRIFUGAL LOAD’ VECTOR 

Refer dgaim to Fie. '2) and .assume ythat tthe .\body ws. rotat- 
ing about the Z axis with aiewlar speed Q.)' Then. the centri- 
fugal force per unit volume at a point distant R from the Z 
axis will be da" R, where p is the density of the material. 


The work done im Eiets "cases 


Ws 027 Ra uediVs . (S50) 
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For constant p22 the corresponding element centrifugal load 
vector is readily obtained by evaluation of the components 
of the integral in Eq. 50 in terms of the nodal variables, 
re. 

il 


1 
eae eS ated GENK.) det oO tieiidedn: (51) 
Je -1 -1 cae ee 


Where det J-is the determinant of the Jacobian coordinate 


transformation matrix (see Appendix Aas 


E. STRUCTURAL BOUNDARY CONDITIONS 
The structural boundary conditions implemented in the 
program are: 
(a2) one or more nodes Prevented from moving axially; 
(b) one or more nodes Prevented from moving radially; 
fe)” the right-hand end Cross-sectional plane remains 
plane and the transmitted axial force is Speciiied. 


Hereinafter this will DewnEeterned tolasmithe plane- 
end boundary condition. 


The computer Program presented in Appendix B has the 
Capability of handling any combination of the above men - 
tioned structural boundary conditions. For boundary condi- 
tions of the types (a) and (b), Simply multiplying the 
corresponding diagonal component of the stiffness matrix by 
1929 gives zero displacement [8] (for practical purposes). 

For the boundary condition of type (c) both ends are initially 
Fixed axially for all solutions. An additional SOlLUEHON 1S 
obtained for unit axial displacement of one end. The axial 
force is evaluated for each solution. Superposition is per- 


formed by adding the displacement vectors for the given 
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loadings (thermal plus mechanical), plus an appropriate 
fraction of the vector found for unit axial displacement. 
immas fraction is chosen so that the resultant axial load has 


mnie specified value. 


Boe S.OLEM EQUATION SOLVER 

Once the desired structural boundary conditions are 
mopbred, then the problem is to find the nodal displacement 
wector 0, corresponding to a given number of load vectors. 


We have 


leas 


OM = SIRS, (52) 


maere K is the system stiffness matrix in banded form and F 
meeceload vector. In the present computer program a single 
Enietesky decomposition is performed on K. Then, by a process 
Se torward and back substitution, each load vector’ is 


meplaced by the corresponding displacement vector. 


feaeeeRINCIPLE OF SUPERPOSITION 

Upon the evaluation of the displacement vectors due to 
the various types of loading, the principle of superposition 
can be applied on the displacement vectors. On each thermal 
displacement vector the displacement due to any other type 
of loading is superimposed and, as the result, the number of 
displacement vectors is reduced to the number of thermal load 


Vectors. 


H. STRESS EVALUATION 
From the system displacement vector 6, the displacement 


vector of each element ce may be obtained easily. Then the 
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eorresponding element strain vector es at any point can be 


found from 


Cae (53) 


1m 
! 
ll 


mnally, the corresponding element stress vector cil is 


obtained by 


eee (ae. (54) 


Its 


Since normally the stresses on the inner and outer surfaces 
moenaxisymmetric bodies are desired, in the computer program 
presented in the Appendix B provision has been made to cal- 
culate the stresses at the two Gauss points corresponding 


Eenc. = + on the inner and outer boundaries of each 


oe 


element (where n = +1). 

Upen the evaluation of the stresses at each point the 
mean stress and the octahedral shearing stress [7] are 
Gameulated. The program gives as output the extreme values 
GE these stresses, the R and Z coordinates of the corres- 


ponding points, and the times of occurrence. 
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V. ONE-DIMENSIONAL TRANSIENT STRESSES 


In this section a one-dimensional comparison of stresses 
is made between exact and finite element results. The tran- 
Sauent temperature problem is the one previously described in 
Bection III. 

ity the slabj edges. are, free) to, translate in) the plane of 
the slab, but are prevented from rotating, the exact solu- 
Eron for thermal stress [9] is 


2 ay ae : 
oy, ao i= (Tove la)? F 5:55) 


where ives the Local Gemperature (Eq. 35), and 


Tavg is the average temperature in the slab 


If we choose E = 2, a = 50: dnd. Ve Ol s¢a lt in consas tent 
units), then the maximum stress obtained from the exact solu- 
ELON 1S 


Ci etd 7780 


mmoere OCCUrs at x = 8, t = 73. 
Using the finite element technique with trapezoidal time 
integration and Irons' correction every 10 steps, the maximum 


eeressaiS.tound. tox be 


Oe ae. = Sets T Le 


widett a21so: OCCUrs at x = 8, +t = 73,.as before. 
Rte1sObserved» that the method chosen gives excellent 


results. 
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Vix TEST PROBLEMS 


Program integrity and accuracy have been verified by 
Solving a number of test problems. Since stresses, whose 
€valuation depends upon derivatives of displacements, are 
known to be less accurate than temperatures, comparisons 
Maen exact results are confined to stresses. Individual 


problems are described below. 


meee THICK CYLINDER 
Monsider a thick cylinder with inside radius 30 inches 


aia sOutSide radius 50 inches and the following material 


Properties. 
Modulus of elasticity Be 26.9 xX 10° Psi 
Borsson's ratio y= .28 : = 
Wecerrecvent Ofsthermal expansion a = 7.22 x 10°° Ee 
See = Btu 

Thermal conductivity kc -=" 28). he ane 
Density op = 489, Lbm/ft° 

we 7 Btu 
Specitic heat eyaen ili homo r 


piearbitrary length of 25 aides has been selected for the 
cylinder and it has been subdivided into two different 
SVeient representations aS in Figs. 4 and 5. The plane-end 
boundary condition with zero axial force is used. The 
Stresses for various types of loading are compared with the 


eerresponding exact analytic solutions as described below. 
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Fecoure 4,9) bworkhadial Elements. Representation 
Of whack Cy lander. 


Figure 5. Five Radial Element Representation 
ot Tnmeck Cy lander. 


1. Thermal Loading 

rom sa Linear Varrativon of the temperature T = Z0R 

ebessStresses Obtained by the finite element method for the 

above representations are compared with the exact analytic 
Solution in Fig. 6. . The Tra for thas problem clearly is 

9 


zero and the one obtained by the program was 10. The 


geeuracy Of the other results is chearly satisfactory. 
a. Pressure Loading 
AUR OTM epressuce, Of 1000 psi.acts on the inner 


Surtace Of the cylinder. Again, t 
10 


RZ is zero and the program 


gives 10 In Fig. 7 the other stresses induced by this 
uniform pressure are compared with the exact solutions. 
Here also the accuracy of the results, even with two radial 


elements, is adequate. 
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FIG. 6 
THICK CYL/NDER UNDER THERMAL 
LOADING, T= 20R 
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Fig.7 THICK CYLINDER UNDER INTERNAL PRESSURE, 
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53. Centrifugal Loading 


The speed of rotation has been assumed to be 500 
mevolutions per minute. The stresses obtained by the pro- 
mean are compared with the analytic solution in Fig. 8. In 
this case also the results obtained by the program, even 
With only two radial elements, are very close to the exact 


solution. 


=. HOLLOW SPHERE 

We consider a hollow sphere with the same material 
meeperties as in the thick cylinder with inside spherical 
radius 30 inches and outside 50 inches. The loading is 
thermal with T = 20 "(spherical radius). Symmetry permits 
Beane Only half of the sphere for the computer analysis. 
whe elements representation is given in Fig. 9. Since the 
preopram gives stresses in cylindrical coordinates, these 
have been transformed to the spherical coordinates for 
comparison with the exact solution in Fig. 10. The accuracy 


et the results is noteworthy. 


©. THERMAL STRESSES IN NOZZLE 

This problem concerns thermal stresses near the inter- 
Section of a cylindrical pipe and the spherical vessel. 
mies tt gives the cross-section of the structure which is 
to be analysed. ihe macerial properties are: 


2853 0° «Psi 


Modulus of elasticity E 


Potsson's ratio VS 550 - 


6 


thermal expansion coefficient a = 7.6 x 10 © 1/F° 
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FIG. 8 
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Thermal conductivity k Oyo -btu/t£t. hr. °F 


B40.5 Lbm/ft~ 


Density 9) 


Specific heat € 128 Btu/Lbm°F 


The loading results from the thermal transient in the fluid 
Meontained in the nozzle. This fluid is in contact with the 
Beructure on surface "1" (Fig. 11) and has the entry temper- 
munre Lime Variation aS in Fig. 12. The fluid in the sphere, 
meren 1S in contact with the structure on surface "2" (Fig. 
i), mas the constant temperature 478°F. At t = 0 structure 
iced Uuilitorm temperature of 478°F- and is stress free. 
Peererror Surface of the structure is tiered: The flow 
metecity past surface 1] is 8.5 ft/sec (inward) and 

—rere is no flow past surface 2. The surface heat 

Beranster coefficients are 1393 and 2910 nvany nee see COs £Or. 
menaces “1 and "2" respectively. 

For the structural boundary condition it is assumed 
that the nodes on the left end cross-section of Fig. 13 are 
prevented from any axial displacement. 

The maximum thermal stresses obtained by the program 


occur in element 14 of Fig. 13 as follows: 


Maximum mean stress = 18.81 ksi. at time = 4 seconds; 
Maximum octahedral Shearing stress = 11.5 ksi. at time 
= 18 seconds. 


These results appear to be reasonable, but no suitable com- 


parison solution is available. 
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O 2 CEC. 
VARIATION OF ENTRY FLUID TEMPERATURE 


IN NOZZLE. (SURFACE 1”. ) 
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SEC. 
FLUID PAST SURFACE ‘2 - 


BiG. 2, Fluid Temperature - Time Histories. 
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VII. CONCLUSIONS AND RECOMMENDATIONS 


A. CONCLUSIONS 

A computer program has been developed for the solution 
Sieaxisymmetric transient heat conduction and thermal stress 
problems. The system will accommodate a wide variety of 
geometric arrangements, thermal and structural boundary con- 
ditions, and mechanical loadings. 

Although double-precision arithmetic is employed through- 
out, efficient algorithms for the manipulation and storage 
of large symmetric banded matrices allow in-core solutions 
wienemodest time requirements. 

ihe quadratic isoparametric elements used allow accurate 
representation of curvilinear boundaries and the stress 
field. Examples presented show that a small number of ele- 
ments is generally sufficient to determine stresses with 


good precision. 


B. RECOMMENDATIONS 

Incorporation of several additional features would 
Sipniticantly increase the program capabilities. The follow- 
ing extensions are recommended. 

1. Material thermal and elastic properties have been 
assumed constant within each element. Since such properties 
are generally temperature dependent, provisions should be 
made for oeriadieavine lie anes them during both temperature 


and stress solutions. 
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2. The surface heat transfer coefficient, taken as 
constant in the program, is a function of temperature and 
flow velocity. Provision should be made to include these 
eerects. 

=. the program presently starts every temperature 
solution with constant initial solid temperature. Provision 
for an externally specified initial temperature vector should 
be included. 

A. The large quantity of temperature and stress results 
generated by the program is currently presented as digital 
printout. Graphical output in the form of two-dimensional 
contour plots of temperatures and stresses should be 


provided. 
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APPENDIX A 


APPLICABLE FORMULAS AND EQUATIONS 


PART 1 


(a): Two dimensional 8-noded (parabolic) 
isoparametric shape functions. 


Corner nodes: 


ey eRe tn) (E, + 1, - 2) 


i O 
Mid nodes: 
¥ . Eee 
g5 = a0), N. = Oy (di S eG a Mod 
i al Pat” 
eco Ne > te el =n ) 
where 
EO = 5 5 > No Fil We 
fb): Temperature at a point in terms of the 


nodal temperatures. 


fe); Coondinates at a point in terms of the 
nodal coordinates. 


(d): The Jacobian coordinate transformation 
Maerix’. 
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(e): Element of the finite element capacitance 
matrix. 


e€ — 
Cae OC J. N; Me dv 


1) 
(f): Element of the finite element admittance matrix. 
ee fe Ne cl 
pay ae 


(g): Elemental volume. 


dViea ek det. de dn 


PART **2 
(a): One-dimensional 3-noded (parabolic) iso- 
parametric Shape functions. 
ia 
End nodes Nee a ee Gis Eo? 
: ) 2 
Mid node N. =p lko =e) 
where, again, Eo = & b3 
(b): One-dimensional capacitance and admittance 
matrices. 
4 23, = a 7 -8 1 ne 20210 
e 2 Per ae eo t 
1% 0 jie IS) Ze i aia 5 8 16 Sime Orr 0 
-1 2 4 in +8 7 Otero 
where & is the length of the element. 
fe): Whe element v vector for the one-dimensional case. 


e ih 


my eel 10 0> 


>/ 


ee eee 


Teh Het = 95 


? ‘fj 
‘fats 
i 
{ 
f 
4 
0 
i 
| 
‘ 
in 
. 


; : PA ik 3 


Lemire fo oxinti sif to Inomalt 
iS 1. A ag = 12 


oly Ts saonall 
mr le he ae 


im ov Igdnemuts 


ene 51 Tomiteq ae 


7 


ivi nn 


veh 


Solano b-saatl 


bend baa 


suboy bigy 


Wrsus ,ored@ 


Ps 
i > ° 
i 
‘i > 
x 
i : .° 
¢ 


PART 3 
(a): Strain-displacement relations. 
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where N. are the same isoparametric shape 
HUMCtlOnN AS in Pamt Gly.(a )). 


(b): The B matrix. 
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(c): The element displacement vector. 
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(id): Elasticity matrix D fOr aleisocropic material. 
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APPENDIX € 


USER'S MANUAL 


In this section the procedure for the use of the present 
computer program is described. It is intended that a person 
With minimum familiarity with the details of the finite 
element method and computer programming be able to use this 
program. 

The steps below are in order and the user is advised to 
follow them carefully. 

For finding the thermal stresses in an axisymmetric body 
indeteaxisSymmetric loading, the user has the choice of using 
either che Metric or British system of units?” “Ihe program 
will handle both systems and the necessary conversions are 
made automatically within the program. Home the unats 
used in each system must be consistent and they should be as 
bisted below in Table V. 

hoOmeLor preparation of the data input fora given axi- 
Symmetric geometry with prescribed thermal and structural 
boundary conditions we go through the details with a simple 


example. 


ocep 1: 
Deaw the longitudinal cross-section of the body to scale. 
identity the cylindrical coordinates R and Z, with the 


erigin of the Z axis on the most left-hand point of the 


1.35 


TABLE V 
UNITS FOR INPUT DATA 


Note: an input card specifies whether British or Metric 


data is being used. 


Mei rake Sst. 


Quantity Bratish ino 
Coordinates ate cm 
Time See. Sec. 
tiemperature °F AG 
Velocity £t./Séc. m/sec. 
Axial Pore lbf kg 
Pressure Lbf£/in? kg/cm? 
Density lbm/ft> gm/cm> 
Specific heat Btu/1lbm°F Cal/ke°oC 
Mog. of elasticity Tee kg/mm? 
Goer. Of thermal exp. dy Aes ESE 


thermal conductivity Beusbr ft. or cal seec.em sc 


Heat transfer coef. Btu/lir. £t.-°F €al/sec.em. “°C 
Load vector lbf kg 
Rotational speed Revolutions/min. Revolutions/min. 


“ibe is pound force and 1lbm is pound mass. 
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mody. The entire cross section must lie in the first quad- 


mamt) of the coordinate plane. 


R 
13. 
2 
5. 
° 5. 10. 15 . 
Pugure: 14."eLongrtudinal Cross Section. 
Beep 2: 


Subdivide the cross section into a maximum of 40 quadri- 
Haceral elements. This subdivision is extremely important 
andewe Should provide smaller elements for the parts of the 
Mm@osecsecctL1on where we expect the laxgest stresses or tem- 
perature gradient. If the body is made from several mater- 
ials, elements should be chosen so that each element contains 
only one material. Any two adjacent elements must share one 
complete side of the quadrilateral. Number the elements, 


Beaneiie trom 1, in any arbitrary manner (see Fig. 15). 


otep 3: 


Simee ‘cight-noded elements are used in the program, iden- 
crtyethiese nodal points around the boundary of each element. 


Four nodes will be at the corners and the other four will be 


MS; 


at the mid-points of the sides. Number sequentially these 
nodes starting from 1 and increasing in the direction where 
the number of elements is the least. (The numbers thus 
assigned are called global node numbers.) For clarity of 
this step assume a cross-section as in Figure 15 such that 
the maximum number of elements in one direction is less than 
the maximum number of elements in other direction. (In 

Fig. 15 we have maxima of 2 and 3 elements in R and.92 soiree 
tions respectively.) Thus we number the nodes beginning in 
the R direction. See Fig. 16. This method will give the 
minimum band width of the stiffness matrix and will save 


computer execution time. 


Figure 15. Subdivision Inot Elements. 


Step. 4: 
At this point we are ready to prepare the figsse. dara 


Eara, iImput quantities are: 


38 


isd sais eff to einiLaq-b 
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Figure 16. -Element and Node Numbering. 


NET The total number of elements 

NNT Total number of nodes 

NCN otal number of corner nodes 

N@EN Total number of mid-side nodes not lying on the 


Strargent linesjoiningethe corner nodes ~(iumber 
Ome OLE” nodes) 


NMAT Number of different materials 


NPRQOB Number of problems to be solved 


For the present we take NPR@B = 1. 


Prepare a single card that reads 
NET, NNT, NCN, N@FN, NMAT, NPRGB 
machethe fLormat (815). For our example, if all the elements 


are from the same material this data card reads: 


Loa ng pa SR mea 


ie) I GOOO0D0NOLHIN9 
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On this card, or any succeeding one, if an input quantity 
(e.g., N@FN) is zero, the corresponding field may have a 0 


or be left blank, unless otherwise stated. 


Beep 5: 

For each corner node a card should be punched which 
meas : The node number I, the R coordinate and the Z coor- 
Gamate Of the node. The format is (15,°2F10.5). 

ies not necessary to sort these cards in the-order of 
increasing or decreasing corner node numbers, they can be 
Pum togcecther in any arbitrary order. <A typical card is 
Menestrated In Step 7. There must be as many punched cards 


as the number of corner nodes (NCN). 


Sscep, 6: 

in this. step. we read in the connectivity array, di w@saye 
for each element we prepare one card which gives the global 
node numbers and the material identification number. 

For each element, start from the lower left-hand node 
and move in the counterclockwise direction within that 
efenene and punch the global node numbers in order. The 
rormat 1S (915) where the last [5 is the material identifi- 
cation number. 

hes very important that the cards prepared for this 
Secpabe put together in order of elements, 1.€., the first 
Gard for element 1, the second card for element 2, etc. 

For ease of sorting the connectivity cards the element 
number may be punched after column 50. As an example, for 


element 2 of Fig. 16 the connectivity card should read: 
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where 1 identifies the type of material in this element and 


Zein, column 55 stands for the element number. 


Seep /: 

If each side of every element is <Piratene. NGFN = O and 
this step is Bar eteds 

For each mid-side node that is on a curved element edge 
wecasa tS prepared with format (15,2F10.5) to read the global 
moae number and the R and Z coordinates. Here, as in step 5, 
mmiese cards may be put together in any arbitrary order. For 
Biemease Of Fig.16, there will be only one card to be punched 


anGedit would read: 


eee tar i Fate 


Q000000000000 o000000 OO000000NNSEHNDNKKODNG EHS eoeeddgaDa 
P2945 67 8 91015 121314 15 16 17 18 19 20.21 22 23.24 25 26 27 26 29 39 31 32 33 34 35 36 37 32 39 49 4) 42 43.44 45 45 47 4B 49 50 Si 52.59 54 $6 56-57 58 
TR te NA eS ie ee ee Iie i DLC fay a Lec | 0 4'a) Asai F ik: 


boa! 
9$ 60 61 £ 


EPS RS SS TT Te ie ies ta [les bite es LK | 


miovassure in all cases,that the double precision capability 
of the program is not compromised by less precise floating 
Pome inputs, it is recommended that all such inputs be made 
in D-FORMAT. The F-FORMAT instructions merely allocate 


input scard fields. 
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ptep 8: 


The properties of the different materials are specified 
teens Step.. For e€ach material a card must be prepared 
euategaves the modulus of elasticity (E), coefficient of 
mieamal expansion (AL), Poisson's ratio (P@I), thermal con- 
dmetivity (TK), density (DENS) and specific heat (SHT). 


fen, read: 
Epa Pol. tke DENS, SHT 


omit ne format, (ZDI2.4, 478.3). The first card will be for 
Macerdal number 1, the second card for material number 2, 


Etc. 


step 9: 

Por this Step a single card must be punched, Starting in 
column one, which reads the word BRITISH or METRIC in accor- 
dance with the system of units used. 

ines completes the input of geometric and material prop- 


erty—data. 


peep 0: 

FOr a transient temperature problem only (no stress cal- 
eumactons) teave a blank card for this step and proceed 
Penecriy £O Step 16. 

FOLea stress problem or thermal stress problem, in this 
Steprwe specity the type of problem and the structural 
boundary conditions. 

The following quantities, when pertinent, are to be 


Specitied. 
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PMEGA Tbeospecad oO: rotation. about the Z axis, in 
revolutions per minute 


PRES The external pressure applied to a boundary 
segment 
TEXT The number of known temperature vectors for 


which evaluation of thermal stresses is desired 


L@GAD ihiesindneatzom ery an. adeitional known load 
vector. If there as one, LOAD.= 1, otherwise 
LQAD = 0 

NNLT Total number of nodes fixed in the longitudinal 
direction 

NNRT Total number of nodes fixed in the radial 
direction 


HPEANE An indication for the plane-end boundary condi- 
fom. ; Lt desiaréed IPLANE = 1 


if any of these items is not applicable, the correspond- 
Mimeeteld 15 left. blank. 
Now we have to punch a single card for this. step which 


reads: 
@MEGA, PRES, ILEXT, L@AD, NNLT, NNRT, IPLANE 


iiehieerormat (2F10.4, S15). 


peep Il: 

it ethere is:°no pressure loading (PRES=0), omit this 
Secu FOL non-zero pressure, the total number of nodes on 
pressure loading boundary segment (NPNT) and the global 
node numbers of these (pressure) nodes (NPN(I)) must be 
Specitied here. Only one segment 1s permitted. 


We read in: 
NPNT, (NPN(I),1=1,NPNT) 


with the format (1215). 
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im-our example; Fig. 16, if there is pressure inside the 


medy, then the data card for. this step will be: 
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y otep i Ae 


iesenere 15 no additional load vector (LQ@AD ="0), omit 
enas Step. 

For LOAD = 1 read in the additional load vector with the 
format (6F10.4). The components of the load vector are 
aemoueed in the order: R component at node 1, Z component 
aeemode 1, R component at node 2, etc. I.e., READ 
(F(1),I=1,NDF), where NDF (the number of degrees of freedom) 


is equal to two times the total number of nodes (2*NNT). 


Sb) ead Gee 

ingany Stress problem, or thermal stress. problem, there 
must be at least one node constrained against longitudinal 
motton. 1.¢.,, NNLI > 0, so.we must specify here which nodes 
are to be constrained against longitudinal motion. These 
global node numbers NNL(I) are read in with the format 


CEES), 1s, ; 
READ: (NN(I),I=1,NNLT) 


In our example, if global nodes 1 and 21 are fixed in the Z 


direction, we have a card that reads: 
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otep 14: 

fe there are no nodes fixed in the ‘radial direction 
(NNRT = 0), omit this step. 

For NNRT >0 we must read in the global numbers (NNR(I)) 
of those nodes which are to be fixed against radial motion. 


iiem.ormat 1s (1015), i-e., 
READ: (NNR(I),1=1,NNRT) 


ieonode 1 and 21 of our example are also fixed in the R 
direction, then the card for this step would read exactly 


aeeee one im Step 13. 


peep 15: 
If there is no plane-end boundary condition (IPLANE = 0), 


Omer this step. 

For the case of end-planes-remain-plane boundary condi- 
elon ise. “(LPLANE = 1), we have to specify the total number 
of nodes of the right-hand end (NNRE) and the global node 
numbers of this end (NNR(I)).- 


We read in: 
NNRE, (NNR(I) ,1=1,NNRE) 


Wien the format (1015). 
Pimetiencase Of Fig. LO, 1f Lt as desired to have end- 
planes-remain-plane boundary condition then the data for 


Eis step will be: 
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Bsani eaiatian aah as 
P000000009 DODOONCOOKODDKONOHNHDOTHNNNHDDNNDOHEHTRIODOHCHHHAOLE 
12345 67 8 $9 1011 1213 1415 16 17 1819 2021 22 23 24 25 28 27 28 29 30 31 32 33 34 35 36 379 hl ede! eg mae hea yen 

eeep 16: 


im this step we begin the specification of thermal 
boundary conditions. Additional details are prescribed in 
meeps 19 through 24. 

Thermal boundary conditions are imposed along discrete 
segments of the boundary. Each such segment must begin and 
end at a corner node of a boundary element. We consider 
Separately the imposition of the various kinds of thermal 
boundary conditions. 

(a) Convection 

Simee the fluid temperature for convection is deter- 
mined from a prescribed temperature-time history at entry 
Section and a specified flow velocity (see Appendix D, 
ZoGckrOn 5), Lt 1S, mecessary that the terms “inside, 1.e., 
adjacent to the symmetry axis, and "outside" have their 
Biarmmary meanings... Thus on the inside surface: the: convec- 
tion boundary condition may be applied to a single (con- 
tinuous) segment. A similar prescription may be employed 
for the outside portion. The entry section used for flow 
ealculations is at the upstream end of corresponding segment. 

(b) Constant temperature 

morspecify constant temperature on a portion of the 
boundary,°the designator “inside™ or “outside” may be used. 


Such a constant temperature portion may consist of several 
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discrete segments. If two different constant temperature 
portions are prescribed, one may be designated “inside” and 
exe other “outside.” 


(c) Combinations of convection and 
constant temperature 


The convection and constant temperature conditions 
Samebe used together, but the portion designated "inside" 
must have only one of these conditions prescribed and the 
Same restriction applies to the "outside." 

(d) Insulated 

All portions of the boundary not included in the 
feements Specifically identified as “inside” and “outside" 
are considered insulated. 

the sto! lowing items, when pertinent, are to be 


paven 2s specified below: 


TINIT The constant initial body temperature 
Q2=0. For insulated boundary condition outside 


o2>0% For convection or constant temperature boundary 
condition outside 


Q3=0. For insulated boundary condition inside 


OS 2.0 a) skOs CONVECTION, Or constant temperature 
boundary condition inside 


OAscOe it Solvine stress problem ‘only 


O4>0:. tf solving temperature problem only 


Q4=0. If solving thermal stress problem 

PPeuvgne whine axial —foree, in. the Z; direction... As usual; 
Pe(plws) for tension and — (minus) for compres- 
Son . 


So, we read: 


MINIES O02, 03, 04, AXIALF 
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mech tne format (SF10.4). 

Paeecase of Fig. 16, if the initial solid temperature is 
60° and we have insulated outside and convection boundary 
Condition inside with 120 kg axial compressive force, then 


Ene card for this step will be: 


be 
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88 


Where the blanks left in the columns. 31 through 40 indicate 


that we want to solve a thermal stress problem. 


step. 17: 

No detion 1s taken in this step - we merely choose between 
preccedine to step 18 or jumping to, Step 26. 

If there are no known temperature vectors for which 
eyaluatiton of thermal stresses is desired, proceed to 
step 18. 

For IEXT>0, i.e., when some known temperature vectors 
ure to beventered for thermal stress analysis alone or 
combined with some other loadings, proceed.directly to 


Seep 26. 


Step 18: 


If (Q4<0), i.e., we are to solve only a stress problem, 


proceed directly to step 27. 
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peep 19: 
If there is an insulated boundary condition outside 
(o25— 0), omit the following steps and start from step 22. 
For the case of a convection or constant temperature 
poungary condition outside (Q2>0.), we have to specify 
the outside heat transfer coefficient (HTC@) (for constant 
temperature HTC@ = ele the Constante, oucside initial 
fluid temperature (TEMP®@) which may be equal to the solid's 
initial cempenmature. the total number of modes im contact 
With fluid outside (NCF@T), and the number of temperature 


ramps for outside flow (NRAMP®). We read: 
HTC@, TEMP@, NCFOT, NRAMP@ 


Weeieehe tormat (D16.4,F10.4,215). ‘See example in ‘step 22. 


pee 20): 

Heneawe. Leadein the global node numbers of the nodes 
icomtact with outside fluid (NCF@(1I)). The sequence of 
these node numbers must be in the direction of flow velocity. 


Read: 
GNGEG Gh), f= ,NCEOT) 


With the format (1215) (see example in step 23). 


Seep ZI: 


Bor cach ramp Of Outside flow we read im the time when 
Ehemnanp starts BDRYO(1,1), the imitrzal temperature of the 
Pano DRYOC(I,2), the final temperature of the ramp (specify 


Meet te ditters from the anitial temperature’ of the next 


ramp) BDRY#(1,3) and the velocity of the fluid for this 


ramp BDRY#(1I,4), with the format (4F10.4) so we read: 
((BDRY@(1I,J) ,J=1,4),1=1,NRAMPQ) . 


The number of cards prepared at this step is equal to the 


“number of Ramps. “see example an step 24. 


Step ey ae 


he ethere iS an insulated boundary condition inside 
(03-0), proceed directly to step 25. 

Ror the case of convection or constant temperature 
boundary condition inside, Q3>0. We have to specify the 
inside heat transfer coefficient (HTCI) (for constant 
inside temperature Hnel=10-") the constant inside initial 
fluid temperature (TEMPI), which may be equal to the solid's 
initial temperature, total number of nodes in contact with 
fluid inside (NCFIT) and the number of ramps of the inside 


flow (NRAMPI). We read: 
HRCr, TEMP? , NCFIT, ~NRAMP I 


miemeche format (D16.4, F10.4,215). 
For our example assume the time variation of entry temper- 
ature £or the inside flow to be as in Fig. 17. The heat 
Peaster Goefticient is 175.0 with initial inside fluid 
temperature 70:0. The card for this: step is: 

a omg i se DU z = 
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Figure 17. Entry Temperature-Time Variation. 


Further details of this temperature-time history are speci- 


mucanin Step 24. 


step 23: 


Here we read in the global numbers of those nodes which 
anemime contact with inside fluid, .(NCFI(1)) with the format 
Cizt>)e, Lt is: important to read these node numbers in the 
orceetton Of £low velocity. 

rer ixample, if the mmside flow :cftFips 116 is from left 


Poeieeit. then £or this step the card reads: 
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Step 24: 


In this step data for inside flow are given. -For each 
famp I we read: the time when the ramp starts BDRYI(I,1), 
the initial temperature of ramp I BDRYI(1I,2), the final 
Pemperature of ramp I (specify only if it differs from the 
freral temperature of the next ramp) BDRYI(1,3), and the 


meer ty Of the,fluid,for this, ramp. BDRY1(1,4).,.We-read.in: 
((BDRYI (I,J) ,J=1,4) ,1=1,NRAMPT) 


Ween tormat (4F10.4). 
Horeexanpie Of Fig. 17, let the flow velocity for the 
first ramp be 15 and that’ for the second and third segment 


be 20. Three cards are needed: 


fan.Dg “Bo0TDo S00.D0 en. ne 
Ed 
3 CARD 
0000 00 0000 ae 00000 00 000002 00 vocgognagsrannn 
1234567 8 $ 10111213 16 | 324 25 25 27 28 29 36 31 32 35 34.38 36 37 3B 39 40 4) 42 $45 6647494 3 t 
100.D0 165.00 eng, pa 2n.be 
nd 
2 WCARD 
000 00 00000000 0009000 00 00002 86 O00000dESSE podeooogs 
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OTH Bo.Da i5.D0 
ofa 
ST CARD 
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W234 6 78) 8 102 131 USMS ABIL IBIAS /202)"22.25'28 25 26 27 28 29 30 31 323 Sa de PTE EGS 4G 47 Jed 51 5253 54 $$ 56 $7 Seale am 
Step ase 


In this step we must specify the following items: 


DTT The time increment (step size) of trapezoidal 
integration 


TIME] The time’ when’ the’ farst° calculated temperature 
vector is to be stored for thermal stress 
evaluation 
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EVERY The interval of time between the successive 
tenperature vectors which are to be stored 


IVEC Total number of temperature vectors to be 
stored for thermal stress evaluation 


INTP tne andicatvon for printing the calculated 
temperature yectors 


ee Neve? 


0 there will be no temperature, primts 


" 


Pe “INTP ] the provcram wai print the tem- 
perature atter every step of tame antegration 


ie INTE) = 5 the program wild print the tem- 
perature agiter every 5S steps of tame inte- 
Pracvons., ELC. 


Mememes Step we prepare a single card that reads: 
Pi wi MiEd a RVERY.: IVEC a NTP 


with format (3F10.4,215). 

foreecxample, af the step size of trapezoidal time inte- 
gration is 2 sec. and we desire to evaluate thermal stresses 
for 10 different temperature vectors, each 40 seconds apart, 
aiemstarting with the first thermal stress calculation at 
meme qual to 60, then the card for this step is: 


moet £4.10 49.00 10 eet 


000000 00 090000 00 OON0OST OOHDNNNNKHHAHHHOTICHOOOAOeE 
8 $10 11 1213 14 15.16 17 18 19 20 21 22 22 24 25 26 27 28 29 20 31 3233 34 35 36 37 3B 39 40 4 42 43 Ad 45 46 47 4B 49 5051 $2 53 54 55 96 57 58 59 6061 6 
Lend 1 1 1 
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Wheres in columns 34 and 35 indicates that the temperatures 
Weltebe printed every 15 steps of time integration. The 
total time of integration for this example is 420 seconds, 


since 


TiME y+ EVERY * (EVEC=1)) = 9420.. 
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Step 26: 

If there are no known temperature vectors (IEXT = 0) 
for which thermal stress evaluation is desired, omit this 
Step. 

For IEXT>0O we must specify the nodal temperatures 
which are to be stored for thermal stress evaluation. We 


read: 
CSHOR BGls Joy <s4<0 =1,,NNT); f=1,, EXD) 


with the format (6F10.4). 

This means that we enter all the components of the first 
nodal temperature vector, followed by the components of 
the second nodal temperature vector and so on until IEXT 


such vectors have been entered. 


Step DT: 


The data cards for this problem are completed now. If 
no more problems are to be solved for the same geometry, 
emat this step. 

If another problem is to be solved for the same geometry 
and spatial discretization, increment by 1 the NPR@B in step 


4 and start the input of the new problem from step 10. 


Step 28: 


i’ there is no access to the IBM 360 computer of the 
Naval Postgraduate School, omit the following steps and 
Sean trom Step 32. 

For the convenience of the user, the author has put 


a so-called CHECK program and the program presented in 
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Appendix B in the computer system at N.P.S. (Naval Post- 
graduate School). It is advised, however, that the user 
check his input data deck with the CHECK program before 
attempting to solve the problem. 

For using the CHECK program prepare the following 


Eonero!l cards. 


//XXXX0000 JOB (0000,0000FT,XX00),'NAME' , TIME=1 

//J@BLIB DD DSN=F0609.BAKH,DISP=SHR,UNIT=2314 , V@L=SER=DUFFY 
//G®. EXEC PGM=CHECK , REGION=100K 

//FTO6F001 DD SYSPUT=A, DCB=(RECFM=FBA,LRECL=133, BLKSIZE=3325), 
// UNIT=SYSQUT 

//FTOSFOO1 DD # 


MWmere the first card is the regular FORTRAN job card used 
Semenrs InStitution. 
Beepare your deck as Pig. 18 and it may be read in from 


the hot card reader. 


| /* 


Data deck 


| control cards of step 28 
| vot card 


EEpure: Los .-Sc€t-up tor Using «CHECK Prop ram. 
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If the output of the CHECK program has any error messages 
or has not been run, then there are some mistakes in the 
Gara Geck or control cards. 

If there are no error messages on the output of CHECK 
program, then study carefully the output and make sure it 
is the same problem that has been intended to be solved. 
Once the correctness of the input data has been insured, 


pmoceced to the next step. 


etep: 29: 


AXITTS stands for Axisymmetric Transient Thermal Stress 
and this is the name given to the program presented in 
Appendix B when stored in the computer. For using that we 


musesprepare the following» control cards. 


//XXXX0000 JOB (0000,000FT,XX00) 'NAME' , TIME=10 

//J@BLIB DD DSN=F0609.BAKH,DISP=SHR,UNIT=2314,V@L=SER=DUFFY 
//GG EXEC PGM=AXITTS,REGION=425K 

//FTO6FO001 DD SYS@UT=A, DCB=(RECFM=FBA,LRECL=133,BLKSIZE=3325), 
a UNET=S¥SGUT ,SPACE=(CYL, (6,1) ) 

PP ETOSEOO1 DD. * 


where the first card is the regular F@RTRAN job card. It is 
advised to ask for 10 minutes time, i.e., TIME=10, since 
ears would not affect the priority of the job within the 
class K jobs.— 

Brepanre the deck aS 1n Fig. 19 (at may be read into the 


eomputer from the hot card reader) and submit a so-called 
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Service request card; available in the computer cénter, to 


ene operator on duty. 


/* 


Data “deck 


Comerol cards Of step 29 


Regular F@RTRAN job card (TIME=10) 


Puowre 9. Neck Set-up for Using AXITTS, Program. 


seep 30: 
Eeaehe user wishes to have a listing of the program he 
May prepare the following control cards as in Fig. 20 since 


Eecceprogram is also listed on the data cell for this purpose. 


Regular FORTRAN job card 
77PRNL EXEC PGM=ITEBPTCH 
//SYSPRINT DD DUMMY 
fy Svsura DD DSN=F0609.BAKH,DISP=OLD,UNIT=2321 , 
Lt V@L=SER=CEL003 , DCB= (RECFM=FB, LRECL=80, BLKSI ZE=2000) 
P7oyxSut2 DD SYS®PUT=A,SPACE=(CYL,6) 
//SYSIN DD * | 
PRINT TYP@RG=PG,MAXFLDS=1,MAXNAME=1 
MEMBER NAME=AXITTS 
REC@RD FIELD=(80) 
/t 


Preuxe 20. Deck Set-up for Obtaining a Listing 
of the Program AXITTS. 
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mie set-up, deck.of Fig. 20 may be read. in. fromthe hot, card 


meader oO get a listing of the program. 


peep 31: 
If the user wishes to obtain a deck of the program AXITTS 


Neemay prepare a deck as in Fig. 21 and read it in from the 
mimeeard reader... He must notify the operator on duty to be 


Prepared for punching almost two boxes of IBM cards. 


Regular F@RTRAN job card 
/ /PUNCH EXEC PGM=IEBPTCH 
//SYSPRINT DD DUMMY 
//SYSUT1 DD DSN=F-0909.BAKH, DISP=$LD,UNIT=2321, 
iM V@L=SER=CEL003, DCB= (RECFM=FB, LRECL=80, BLKSI ZE=2000) 
feorsut2 DD SYSOUT=B,SPACE=(CYL,6) 
fpovoIN DD * 
PUNCH TYP@RG=P%,MAXFLDS=1 ,MAXNAME=1 
MEMBER NAME=AXITTS 


REC@RD FIELD=(80) 


/% 
Pogeure 21. Deck Set-up for Obtaining a Punched 
Deck of the Program M1 Tis 
Seem OZ: 


Foreusane the program, presented an Appendix B, if the 
Mserraoes wot. have access to the computer facility at N.P.S., 


he must punch a copy of the program. (Good luck!) 
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For the storage location 425K bytes are required and 
Since the output may be longer than what usually is allowed, 
three additional ev intens are recommended. The execution 
time required depends on the size of the problem and the 
nomper ‘of time integration steps; however, 10| minutes 
computer time would be sufficient for a problem of 140 nodes 
and 500 steps of time integration, with 20 temperature 


weecrors stored for stress calculation. 
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APPENDIX D 


PROGRAMMING 


The computer program presented in Appendix B is formed 
from one main program and fifteen different subroutines. 
The communication between the main program and the sub- 
Houcines is handled through the common blocks. In order to 
minimize storage requirements, most of the storage location 
Gietme system stiffness matrix is used for temperature 
Gareulations. This is accomplished by use of EQUIVALENCE 
ents. ihe total storage: requirement is 425K bytes. 

mimes section the funetion Of Some parts of the pro- 
gram is discussed and the assumptions used are brought to 


attention. 


1. MAIN PROGRAM 
The main program is simply calling different subroutines 
when they are needed. The subroutines which are called in 


main program are: 


INPUT, PR®@B, CANDY, FL@W, FORMV, TEMPER, PL@T, STIFF, FO@RMF, 


CENTP, PRESS, DISPL, AND STRESS. 


72 SUBROUTINE INPUT 

In this subroutine the data regarding the geometry of 
the body, subdivisions into elements, and the material 
properties are-read in and printed ie Catenlataon Of che 


Mid-side coordinates based on the straight line is done 


o t1omasdh — 
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' in f Cc #2e7an(bToeG? abt 
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mere. This subroutine calls for Subroutine PLOTRZ once 

to produce a graphical representation of the nodal points 
then the half-band-width of the system stiffness matrix is 
calculated from the connectivity array. Also, depending 
upon the choice of the system of units, the necessary con- 
versions in each system of units are accomplished and the 
mererence temperature (70°F or 20°C) based on the system of 


units is selected here. 


3. SUBROUTINE PR@B 

Por each problem this. subroutine is called once by the 
main program to read in and print out the information 
regarding the nature of the problem. Based upon the given 
information the type of the problem is distinguished here 
and for the thermal problems the length of time integration 


momcalculated. 


4. SUBROUTINE CANDY 

Subroutine CANDY (C and Y) evaluates the capacitance 
iacrrx C and admittance matrix < Of Eq. 3, an banded form: 
Eereecen problem this subroutane 1s called once. The func- 
pional —tlow chart of the-subroutane GANDY is given in 
Pieces oince the only non-zero elememes of the symmetric 
matrix ‘a (fq. 11) are the, diagonal, t07st and second super-= 
Geagomals. in order to conserve Storage and reduce the num- 
bec Of arithmetic operations, the non-zero elements are 


stored in three different vectors (DIAG, @FF1, @FF2). 
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Fieune ZZ. Functional Flow Chart of CANDY. 
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Once the system capacitance and admittance matrices C 
+ 
wee are formed, then the matrices A and G (Eq. 29) are 


evaluated and replace € and Y respectively. 


5. SUBR®UTINE FLOW 

If there is any convection or constant temperature 
thermal boundary condition, this subroutine is called from 
Eve main program for each step of time integration in order 
bemevaluate the temperature of the fluid nodes. The calcu- 
ieron 1S based upon the constant fluid flux assumption 
forepeth inside and outside flow). Consider a section of 
pimenresular cylindrical pipe as Fig. 23 and focus attention 


on element i+] on the inner boundary. 


Becunres 25. Fluid Node Representatron, for inside Flow. 


fe fluid is moving from left to right we may number the R 
Ene e coordinates of the fluid nodes of the element i1+1 as 
imebao. 25, then the volume of fluid pumped into the pipe 


aeotane £ 1S given by 
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Qin =) arV is (56) 


where a and v are, respectively, the area and velocity at 
ehewentrance section. Also the volume of the pipe up to 
the mid-node of element itl is 
2. : 
Qa iz Va Poh Ti Wee (7) 
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where Vv. is the volume of the pipe between the entrance 
section and element i+1 on the inner boundary. 


Now we can write 
3 
eae) RN. (58) 


where N; are the one dimensional shape functions (given in 
Appendix A). 
fiiewe Substitute Eq. 58 into 57 and make the coordinate 


transformation, after integration we get 


_ 1 : 2 DD : : 
Qy = Vat zyg (22727) (S1R]+64R5+RZ+46R,R,-8RyRz-14R,R,). (59) 


Noweby equating Eq. 59 to Eq. 56 at any time t we may deter- 
mine whether the front of the flow has already passed the 
mid-node of the wetted side of the element itl. 

Mcoinpbar arsument can be carried out for the end node 
Of the element i+1l. In that case we also get an equation 
wma tOubaG. 59 with different numerical coefficients. 

For the case of the outside flow it has also been 


acsumed to have a constant fluid flux and a fictitious 
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entrance circular cross-sectional area (whose radius is the 
largest R plus unity) has been assumed. 

Phe numerical coefficients! of R'srin Eq. 59 are used 
in subroutine FLOW and it has been assumed that the temper- 
ature of a fluid “particle” does not change during passage 
through the active section. This is believed to be an 
acceptable approximation for representative values of flow 
welocity and active length. With this assumption, for a 
Geuven entry time-temperature relation (inside and outside 
flow) this subroutine evaluates the fluid nodal temperatures 


at any time. 


6. SUBR@UTINE F@RMV 

me every step of time Integration this «subroutine is 
eameca by the main program ‘to form the vector v of \the 
peone-hand side of the discretized finite element (Eq. 153) 
for the given thermal boundary conditions. The program 


iecet 15 Self explanatory. 


7. SUBRQ®UTINE TEMPER 

ines nodal temperatures are calculated in this subroutine 
werieche werapezordal time integrations. . Irons’ correction 
aeapplved here after every L0“steps of tame integration: 
fice tenperacure vectors selected) for the stress) analysasare 
Stored. The transient temperatures will be cprainted out at 


the desired interval of time. 


165 


¥ 
_ 


5 = 
1400") Sets LaRG lise e=slGrg Selene) 


Lonaze 169d Sen ty rin” wit a2 

i ai e' Yo etitecsstieoy tageseee 
weir eeod 2ad ae bas WAT sat rt 

a 


io 107 Bo0h “ofskssaa" Braee 


_ 
: (aT .noliae2 syvigos ort 


(oz. Mortem indtjghy eee 


14n8l. svitae bie 


julaToOgmhal-on7 Ae 


sis! > Snituowlie ea 


7 . ees 


: 
AGT (Ta 
lO tore 
| inn 44 
b al Pia 
; ; ne 


8. SUBR®UTINE STIFF 

For any stress or thermal stress problem subroutine 
STIFF is called once by the main program to evaluate the 
stiffness matrix of the system. Once the stiffness matrix 
of the system is calculated then the desired structural 
poundary conditions are applied and, at the end, the Cholesky 
Pecomposition is performed. The functional flow chart of 


Supneutine STIFF is given in Fig. 24. 


9. FQRMF 

In subroutine F@RMF the thermal load vectors are calcu- 
lated for as many as (IVEC) given temperature vectors. 
these thermal load vectors are the columns of a rectangular 
Hiatrix F. The provision is made that no matter how many 
eEemperature vectors are given (always EVEC < 20) the IVEC 
pomber of columns of the F dhe fribtedawich the thermal load 
wMecuers and the last column of F will comtain the load 
WeeEor corresponding to’ the unit end displacement for zero 
evieecoree when the plane-end boundary condition is applied. 


thesttow chart of subroutine FORMF is ’eiven in Fig. 25. 


HOS SUBROUTINE CENTE 

Eee chessystem 15 rotating about the axis of revolution, 
his Subroutine is called once by the main program to eval- 
mace te centrifugal load vector.s'This load vector is 
minds placed in the first columm after the thermal load 


weerors (VEC! position) in F< 


166 


Fa ‘A 
' y ergoyy Reaw eng Nd word fe 


izo6b eft gedit bade tipofes a nateye 


1 bite Bei Tage 


i, 


bemrsbiaq-e8 je tsuee 


nin 


ity SMING enisueveae 


af _ - 
'r e3n0- .nSdeye o> Toes 5M 2m 
wtb enol iig 
f ‘9 


(O9V1) 28 Yann 26 
Af | 


aaa 
A . i 
(due woldovd erage femokn 1% ceothe 


paiva = 
has 
heey: 


“F 


7! 
2) 
-_ 


et 1416 ame 
¢ 


rorsove wel brea 
: icntvo rg pili ee 
| hf £40 ™ ‘Vv " 7 

ii to ai Poa 
> amt dae had 


und bow é ogee 


lauvria2f6> ae 


dv sate 


‘to S7ang 
< 


wo ey TUG 
wuge vr 


ang 


aN (tues | 


i bene 


ent | 


P+O040) 


L=1,NET 


Form decomposed D matrix 


Element stiffness 
matrix += 70 
I=1,NGP 


GX (1 
ae 


N = 
6 = GX) 


Evaluation of shape 
Function and 


an aN 
Je "2 son 


Formed 2acobaan 


Form B Ma tx 


Evaluate element 
Staciness matrix 
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| 
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conditions 
Decompose system 
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RETURN 


Rueure 24. Functional Flow Chart of STIFF. 
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Figure 25. Functional Flow Chart of FQ@RMF. 
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11. SUBR@UTINE PRESS 

If there is any pressure acting on the system, the 
Pressure load vector is calculated in this subroutine and 
Enesevector is placed in the column next to centrifugal load 
weetor, if any, otherwise it will fill the position desig- 


maeea tor centrifugal load vector in F. 


2. SUBRGUTINE DISPL 

ihemduspiacement’ veetors are found im this subroutine. 
Since the stiffness matrix is already decomposed in banded 
romm; then ene displacement vectors one after another are 
obtained by back and forward substitution. . The principle 
Gf Superposition is applied here and finally the axial 
Force, if any, 1S corrected for the plane-end boundary 


condition. 


ie SUBROUTINE STRESS 
For every problem this subroutine is called by the main 


pageram once to evaluate the stresses at the points 
1 


v3 
caleulated and printed for each displacement vector. For 


(ia2l, b= 4 Vor "lachvetencnit. Mhesiransient Stresses, are 
each problem the maximum mean stress and the maximum octa- 
Hhedwal shearing stress; together with the corresponding 


tames and locations, are printed. 


14. LIMITATIONS 
In the process of the development of the program.das- 
eussea so far, it has been intended that al of the calcula- 


teens and storage of data occur in the computer without 
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using any external devices such as disks or magnetic tapes 
so as to be able to solve any sizable problem in relatively 
short time. 

pore this reason the following limitations’ (Table Vi) are 
set forth which give an overall size of 450K bytes to the 


program. 


TABLE VI 


MAXIMUM VALUES FOR PROGRAM PARAMETERS 


NBAND Half band-width of the system stiffness matrix 66 


NMAT Total number of different materials 5 
NET Total number of elements 40 
NNT Total number of nodes 149 
NCFGT Total number of nodes in contact with Si. 


outside flurd 


NEPIT iota number Of modes in contact with j Soy) 
imnsade fluid 

NNRE Total number of nodes of right-hand end By 

NNLT Total number of nodes constrained against any Su 


longitudinal motion 


NNRT Total number of nodes constrained against any 37 
radial motion 


VEC Total number of temperature vectors stored for 20 
thermal stresses 


NRAMP® Total number of ramps for outside flow WS 
NRAMPI Total number of ramps for inside flow 15 
NPNT Total number of pressure nodes 57 
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